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SECTION 1 


INTRODUCTION 


The purpose of the work reported here was to develop a computational strategy by 
which spacecraft thermal designers could devise an optimal thermal control system to 
minimize thermal distortion. Control would be maintained by the optical coating 
pattern chosen for the external surfaces and the heaters chosen to supplement the 
coatings. The coatings and heaters were to maintain the required temperatures, 
temperature differences, changes in temperatures and changes in temperature differ- 
ences for specified equipment and elements of the structure. A strategy appropriate to 
computer-aided design was anticipated. 

Most spacecraft use optical coatings and heaters for temperature control. The optical 
coatings are placed on the external surfaces to control the amount of sunlight absorbed 
and the amount of infrared radiation emitted. If the temperature control obtained with 
the coatings is adequate, the design imposes no weight or power penalties on the 
mission. However, for the increasingly sophisticated mission-related equipment, such 
as optical telescopes, the temperatures must be controlled closer than is possible with 
coatings alone. The uncertainty in the coating properties and the changes in the 
coating properties over the mission life usually result in an unacceptable temperature 
range. In addition, the vehicles must usually operate in many different orientations, 
with respect to the sun, so that a coating pattern appropriate for one flight might be 
totally inappropriate for another. Therefore, the coatings must be supplemented with 
another temperature-control method. The method normally chosen is thermostatically 
controlled heaters. These are accurate, light, reliable, and easily applied to any part of 
the structure or equipment. However, heaters draw power and therefore impose both 
weight and power penalties on the mission. 

The pioneering work in the optimal thermal design of spacecraft was done by Costello 
and others in 1968 (Ref. 1). Temperatures were to be controlled by coatings alone. In 
addition, only the time-average temperatures were considered, with the time period 
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being the orbit period. Casagrande attempted to extend this work to inciude the 
fluctuations in temperatures during the orbit (Ref. 2); however, he was successful only 
for ideal cylindrical satellites with one internal node. Neither investigator considered 
the uncertainties in the coating properties; therefore, neither method has found much 
use in industry. 

In the present work, we focused on the development of a practical tool for the designer: 
an optimization program compatible with his present analytical tools and useful in the 
conceptual-design phase. We initially examined the general problem of spacecraft 
temperature control, including all phases of the mission and all phases of the design 
process. Indeed, the optimization criteria developed are generally applicable. 
However, the primary need appears to be in the conceptual-design phase, when the 
weight and power penalties have an important impact on mission planning. After the 
weight and power budgets are established, optimization is less important; it usually 
consists of maximizing the design margins within the budgets. 

We have succeeded in developing an optimization strategy including all of the elements 
normally considered by the designer: various vehicle orientations, various equipment 

power dissipations, coating-property uncertainties, uncertainties in the internal thermal 
conductances, and coating-property changes. 

The method uses the same thermal model data the designer develops for his analyses; 
there is no extra work. He can use his intuition to speed the optimization process and 
to simplify the resultant heater arrangement, but he can also use the process as is. We 
have demonstrated how the optimization process can be computerized. The results will 
be in the form of a specification of the coating pattern to be used on each external 
surface, the heater capacities required for each node, the average power required for 
each mission, and the maximum and minimum temperatures expected for each node. 
However, the method is limited to quasi-steady-state temperatures. The designer must 
select time periods over which steady-state temperatures are representative of the 
mission temperatures. Because most sensitive equipment is located well inside the 
vehicle, his selection is usually not difficult and the orbital period is adequate. 
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SECTION 2 


BACKGROUND 


The complexity of the mathematics describing the optimization problem presented in 
the following section, gives some indication of the complexity of the problem facing the 
thermal designer. There are too many variables to permit a manual, systematic search 
for the optimal design. Usually the designer will start with a simple model of the 
vehicle with a uniform external coating pattern, and compute the resultant 
temperatures. He may vary the coating pattern to offset some of the extreme 
temperatures (if he can foresee the effects for all possible missions and vehicle 
orientations) or he may select a uniform coating pattern that keeps the equilibrium 
temperatures low, and then select heaters to maintain the desired node temperatures. 
As the vehicle design becomes more definite, he builds a more detailed thermal model 
of the vehicle, computes the temperatures for representative missions, and modifies the 
design — usually as little as possible — in the search for a coating pattern that keeps 
the heater power within the budget predicted with his original simple model. Most of 
the design modifications are in the form of interior insulation, thermal straps, special 
component mounting procedures, and, if necessary, heat pipes. Usually it is too 
difficult to pick an external coating pattern because the thermally most severe missions 
or vehicle orientations then become difficult to identify. 

The Coating Selection Program developed by Costello under NASA contract in 1968 was 
an early attempt to assist the designer (Ref. 1). Initially the program was used 
extensively. Some of the early optimizations indicated that some spacecraft for which 
the designer selected a black coating would have performed better thermally with a 
polished surface. However, for tight temperature control, the uncertainty in the 
properties of the polished surface gave unacceptable uncertainties in the equipment 
temperatures. The Coating Selection Program soon fell out of use and the designers 
went back to their semi-intuitive methods. 
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To avoid the shortcomings of the Coating Selection Program, we interviewed many 
thermal designers and some systems engineers. Our own experience in thermal design 
helped. From the interviews, we developed a list of elements that should be considered 
in the thermal design process. From these, we synthesized a strategy that included all 
of the elements. 

The two-tier aspect of the design process was an important consideration in developing 
the strategy. Early in the design process, the weight and power budgets are undefined. 
The thermal models usually consist of at most several hundred nodes, and design 
tradeoffs can have an important impact on the mission. Later in the design process, the 
designers seek to increase the design margins within the fixed budgets. We have chosen 
to develop the thermal-optimization procedure for the early stage of design because at 
that point it has a greater impact on the mission. 
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SECTION 3 


DEFINITION OF THE OPTIMIZATION PROBLEM 

3.1 SIMPLE EXAMPLE 

Before presenting the precise mathematical formulation of the optimization problem, 
we will present a sample problem. The rigorous mathematical definition of the 
optimization problem tends to become confusing because so many symbols are required. 
The sample problem will help in understanding the more abstract equations presented in 
the next section. 

Consider the three-node problem shown schematically in Figure 1. The vehicle is 
divided into two external nodes, Node 1 and Node 2, and one internal node, Node 3. The 
vehicle can be envisioned as two concentric cylinders that are infinitely long. 



Figure 1. Three-Node Problem 



To keep the problem simple, we will linearize the radiation terms. The heat-balance 
equations for the nodes can then be written: 


0 = K 13 • (T3-T l) + K12 ' (T2-T l) + al * Si + ej • (E1-H1 • Ti) + Qi 

0 = K 2 3 • (T3-T2) + K12 • (T i-T 2) + a2 * S2 + e2 ♦ (E2-H2 * T 2 ) + Q2 (1) 

0 = K13 • (T1-T3) + K23 * (T2-T3) + Q3 


Where 

Kjj = the conductance between Node i and Node j 

Tj = the temperature of Node i 

aj = the solar absorptance of Node i 

Sj = the solar radiation on Node i 

ej = the infrared emittance of Node i 

Ej = the earth radiation onto Node i 

Hj = the thermal conductance to space from Node i 

Qi = the electrical energy input to Node i either as heater power or as 

equipment dissipation. 

The a\ and e[ are determined by the patchwork of coatings put on these nodes. For 
example, if a solar absorptance of 0.6 is desired, 50% of the surface might be coated 
with black paint and 50% with polished aluminum. The corresponding emittance might 
be 0.55. In what follows, we strive to find the percentages, or fractions (fj) of each 
surface that should be covered by each of the several coatings available for use on 
spacecraft. These fractions are the independent variables, or design variables, of the 
optimization problem. They will be chosen to minimize the weight associated with the 
heater power. 

3.2 GENERALIZED PROBLEM DEFINITION 
3.2.1 Objective Function 

We want to minimize the weight penalty (W) of the thermal control system, which is the 
weighted sum of the heater energy required for the most severe environment (Q), and 
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the weighted sum of the heater capacities (C n 's). The weighting factor for Q will be 
the weight of the power system, including the solar cells and batteries, as required to 
meet the maximum energy demand. The weighting factor for the heater capacities will 
be the weights of the heaters and heater controllers as well as the weight of the 
harnessing required to carry the power to the heaters. The C n *s could depend on n, the 
node number, because the length of the harness could depend on n. In equation form, 
the function to be minimized is: 


W = j^max ab • Q + bb (C n ) + ae • Q + be • £(C n )J 

where 


( 2 ) 


Q = max [£<q nm >] = maximum over all missions, m, of the total power 

demand, summed over all of the nodes, n 

Q nm = orbital-average power demand for node n in mission m 


C n = max (Qnm) = heater capacity for node (n) which is equal to the 
maximum power demand over all missions (m) 

ab = weight that must be added to the power system to provide an additional 
unit of the power, based on power demand at the beginning of the 
mission 


ae = weight that must be added to the power system to provide an additional 
unit of power, based on the power demand at the end of the mission 

bb = weight of the heaters and associated controls, which is usually propor- 
tional to the heater capacity, based on the power demand at the 
beginning of the mission 

be = weight of the heaters and associated controls, which is usually propor- 
tional to the heater capacity, based on the power demand at the end of 
the mission. 


The weight penalty must be evaluated based on conditions at the beginning and at the 
end of the mission and the maximum penalty selected as the optimization criterion. At 
the beginning of the mission the output of the power supply, such as solar cells, is high 
because there have been no failures and no deterioration of the cells or batteries. For 
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example, at the beginning of the mission, the power system usually supplies 32 volts. 
At the end of the mission the voltage can be as low as 22 volts. Therefore, a high 
demand for power at the end of the mission is more serious than at the beginning. 
However, the heaters and controllers must be sized to handle the maximum voltage. 
Therefore, both phases of the mission must be included. 

There is usually a fixed amount of power and weight associated with simply having a 
heater, independent of the size of the heater. Some power is required to sense the 
temperature (independent of the amount of power supplied). The weights of the sensor 
and the controller are nearly independent of the size. However, if these fixed weight 
penalties were added to W, they would not vary with coating patterns and heater 
locations; therefore, they would not enter the optimization process. Consequently, 
there is no need to carry these terms in the expression for W. 

3.2.2 Design Variables 

The design variables in the optimization process are f nc , the fractions of the area of 
node (n) that are covered with coating (c). For each coating pattern, the heater 
capacities and powers are determined by the heat-balance equations and the tempera- 
ture restrictions, all of which are constraints placed on the optimization process. The 
restriction on f nc is 

0.0 < f nc < 1.0 (3) 

The coating properties consist of the solar absorptance (a) and the infrared emittance 
(e). These are subject to a range of uncertainty, (da) and (de), and to changes, (ca) and 
(ce), from the beginning to the end of the mission. Two sets of nominal temperatures 
must be computed for the beginning and the end of each mission. Both sets must be 
tested against the design requirements. In addition, the uncertainties will result in 
uncertainties in the node temperatures; therefore, each computed nominal temperature 
must be re-computed to include its range of uncertainty. The upper and lower limits of 
this range are the temperatures to be used in determining whether the node tempera- 
tures will be held within the allowable limits. 
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3.2.3 Constraints 


The design constraints are 

o The temperature of each node (n) during each mission (m) must be within 

the prescribed lower and upper limits (Ln) and (Un)* The weighted-average 
temperatures of each designer-defined group of nodes must be within the 
prescribed limits (GLn) to (GUp)- 

o The temperature difference between each pair of nodes (n) and (j) must be 
within the prescribed limits (L n j) to (U n j)- 

o The weighted-average temperature differences between each designer- 
defined group of nodes must be within the prescribed limits (GL n j) to 
(GU nj ). 

These conditions are written: 


f_T nm U n 


(4-D 

GLn — Y (WTn * Tnm)/N 

GU n 

(4-2) 

Lnj < T n m - Tj m < U n j 


(4-3) 

GLnj <_Y (WT n * T nm )/N - 

^ (WTj • Tj m )/3 GUnj 

(4-4) 


where WTn is the weighting factor for node n for computing the weighted average 
temperature of a group of nodes. Each of the T's in the foregoing constraints must be 
interpreted as the nominal temperatures for the beginning and the end of the missions, 
plus or minus the uncertainties. 

In addition to these design limits, the First Law of Thermodynamics must be satisfied. 
We can write the heat balance equation for each node (n) in the form 

Y Q^nj * (Tnm-Tj m )J +£ jj^nj * ( B nm - BjmC] 

= Qnm + D n m + Y ( a c * *nc) * ^nm * ^n + ]T ( e c * *nc) (5) 

* An * (Enm'^nm) 
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Where 

B nm = blackbody radiation at temperature Tnm 

T nm - temperature of node n in mission m 

^nj = conductance between nodes n and j 

R nj = radiant interchange factor between nodes n and j 

Q nm = heater power for node n in mission m 

D nm = equipment electrical heat dissipation into node n in mission m 
a c = solar absorptance of coating c 

f nc = fraction of the external area of node n covered by coating c 
S n m = solar flux on node n in mission m 
A n = external surface area of node n 

e c = infrared emittance of coating c 

E nm = infrared environmental flux incident on node n in mission m. 

The sums on the left-hand side of Eq. (5) are over the j nodes that are thermally 
connected to node n. The sums on the right-hand side of this equation are over the c 
coatings that might be applied to node n. 

The heat-balance equation is also used to determine the uncertainties in the node 
temperatures as a result of the uncertainties in the coating properties. If the coating- 
property half-ranges, da and de, correspond to three standard deviations so that they 
encompass 99.6 percent of the possibilities, then the corresponding uncertainties in the 
temperatures can be computed as the root-sum-of-squares (rss) of the uncertainties due 
to each coating-property uncertainty. Therefore T in the foregoing constraint 
inequalities must be interpreted as the temperature, determined from the nominal 
coating properties, plus or minus the uncertainty in the temperature. The more 
stringent of the plus or minus limits is to be used in satisfying the inequalities. 

With a manageable increase in computational complexity (but possibly a burdensome 
increase in input requirements) the optimization process can be enhanced to include 
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failed missions. Reduced temperature-control requirements could be included if the 
acceptable temperature limits (Ln, Un, Gn, GUn> e tc.) were made mission dependent. 
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SECTION 4 


COMPUTATIONAL METHODS 


The practicality of the optimization strategy depends on the computational cost. There 
are so many complex computations that a computer is clearly required. However, even 
with the computer, the computations could be too lengthy to be practical. Therefore, 
we have examined the computational methods in use today and typical solution times so 
that we can estimate the cost of running the optimization strategy on a computer. 

Four types of computations are required: 

e Solving the heat-balance equations 

• Determining the effect of uncertainties on the temperatures 

o Computing the derivatives of the objective function with respect to the 
design variables 

• Selecting the next set of design variables on the path toward the optimum 
design. 

Each type is considered in the following paragraphs. 

If there are M missions, N nodes, P models based on different internal conductances, C 
coatings, K heaters/ temperature-sensitive nodes, and I iterations to find the minimum, 
then the number of equations that must be solved is as shown in Table I. To help in 
recognizing the dominant problem, we have included in Table I the number of 
computations required for a typical problem containing 20 missions, 150 nodes, 4 
models, 5 coatings, 30 heaters, 30 temperature-sensitive nodes, and 25 iterations. Note 
that the number of computations doubles because all cases must be considered at the 
beginning and at the end of the mission. It is evident that the uncertainties and the 
derivatives of the uncertainties dominate the computational effort. 
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TABLE I 

THE NUMBER OF SIMULTANEOUS EQUATIONS TO BE SOLVED 


Set of Equations 

Number of 
Equations 

Number of 
Times Solved 

Typical Number 
of Times 

Heat Balance 

N 

2MPI 

4,000 

Uncertainties 

N 

4MPCI 

40,000 

Derivatives of W 

2C 

I 

25 

Derivatives of T 

N 

2MPKI 

120,000 

Derivatives of dT 

N 

2MPKCI 

600,000 
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4.1 SOLVING THE HEAT-BALANCE EQUATIONS 

A review of the available computer programs reveals that the following methods for 
solving simultaneous equations are in common use: 

1. Gaussian elimination 

2. Successive over-relaxation (SOR) 

3. Gauss-Seidel 

4. Conjugate gradient 

5. Choleski factorization 

6. Matrix inversion. 

All operate on the linearized equations and require iteration on the non-linear terms. 
The Newton-Raphson technique is the most popular method for linearizing the 
equations, although the resultant set of equations is no longer symmetric. Of the six 
solution techniques, Gaussian elimination is the most popular. No commercially 
available programs use matrix inversion, although matrix inversion was found to be the 
most economical (Ref. 1). 

The data Costello generated in Reference 1 showed that matrix inversion requires 
approximately one-sixth of the computation time required by SOR. Therefore, its 
running times are comparable to the incomplete Choleski (IC) methods. In addition, 
matrix inversion gives the derivatives that are required for the optimization. Thus, it is 
the most promising of the solution methods. 

The data-storage requirement of matrix inversion had previously been considered 
prohibitive and the accuracy of the inversion is a constant worry. However, with 
present-day computers, especially those with virtual memory, the storage requirements 
are less important and inaccuracies can be circumvented by iteration. 

If we anticipate the problem of determining the uncertainties in the temperatures, as a 
result of the uncertainties in the coating properties, we can eliminate most of the 
methods. To determine the uncertainties, we must calculate the derivative of each 
temperature with respect to each coating property. The uncertainty in the temperature 
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is equal to the root-sum-of-squares of the derivatives, multiplied by the uncertainty in 
the corresponding coating property. In effect, we must determine each term in a 
matrix, square each term individually, and sum the rows. 

The only methods that need be considered are those which yield the derivatives 
directly. Consider a method that does not yield the derivatives directly, such as 
Gaussian-elimination. Counting the solutions for the beginning and the end of the 
mission, the heat-balance equation must be solved 2MPI times (4000 times) to 
determine the temperatures. For each of the 2MPI times, the derivative of K 
temperatures with respect to 2C coating properties (one solar absorptance and one 
infrared emittance per coating) must be calculated so that the uncertainties can be 
determined. Therefore, 2KC (300) derivatives must be computed. If the derivatives are 
determined by Gaussian elimination, the derivatives of all temperatures with respect to 
the 2C (10) coating properties are determined simultaneously. Therefore, 4MPIC 
(40,000) sets of N (150) equations must be solved in the course of a solution to the 
optimization problem. 

For the Gaussian-elimination method to be competitive with the matrix-inversion 
method, it must be 2C times faster than the matrix method. The results in Reference 1 
indicate that the two have comparable speeds; the matrix-inversion method is faster if 
the inverse can be determined only once. Therefore, Gaussian elimination will require 
approximately a factor of 2C longer running time. 

Of the methods listed, only the Choleski and matrix inversion methods can produce the 
inverse matrix. The Choleski method is limited to symmetric matrices; it is probably 
faster than direct matrix inversion. Therefore, the choice of the two methods depends 
on whether we want to restrict the problems to the most common case in which the 
conductances are symmetric. This restriction would prevent the use of the method to 
problems with fluid flow. This restriction would also prevent the use of a one-way 
conductor, as permitted by the most popular public domain thermal-analysis program, 
SINDA. One-way conductors are artificial, so that little would be lost by imposing this 
restriction; however, this restriction may require the user to modify his thermal model 
and would thereby violate one of the desirable characteristics of the optimization 
program: not adding to the modeling burden of the user. The limitation to non-flow 
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problems is sufficiently restrictive to conclude that the matrix-inversion procedure 
should be used. 

To solve the system of equations by matrix inversion, we recast equation (5) in the form 

^ d^nj * (Tnm-Tjm)^j = V nm +£ £pnj ’ (Tnm - Tjm)TJ (6) 

Where 

^nm = Qnm + Dnm + £ fee ' *nc) ’ ^nm • A n + ^ (e c * f nc ) • A n • (E nm -B nm ) 

“Z [j^nj * ^nm'Tjm)^] [j*-nj * (B nm -Bjm)|J . (7) 

Note that Vnm is the residual of Eq. (5). When V nm is zero, T nm is the solution. The 
equation is solved iteratively by assuming T nm , computing V nm , and solving Eq. (6) for 
the new T nm (on the left-hand side). The solution is obtained by inverting the matrix of 
coefficients (G n j) which is equal to the conductance (K n j) plus the linearized radiation 
coupling. According to the study done in Reference 1, the computation time is 
minimized if the linearization is done symmetrically. Then 

Gnj = Knj + Rnj • (Tnm’Tnm + Tjm'Tjm) • (T n m + Tjm) = Knj + RLnjm (8) 

In addition, the same inverse can be used for all missions if the temperatures are chosen 
to be the maximum expected for all missions. Therefore, the inverse needs to be 
computed only once for a given optimization. Updating may be beneficial after several 
sets of design variables are evaluated. An approximate inverse might also be useful. 
The approximate G n j for mission (m) is 

Gnjm = K n j + RL n j + RL n j m - RL n j (9) 

Where 

RLnj = the second term on the right hand side of equation (8) with the 
temperatures equal to the maximum 

RL n j m = the second term on the right hand side of equation (8) with the 
temperatures equal to those of mission m. 
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We can write 


Gnjm = (Knj + RLnj) * ^nj + (Knj + RLnj) - ^ ' (RJ-njm - RLnj)^ (10) 

where 

I n j = the identity matrix. 

If the last term is small, we can form the approximate inverse of G n j m 

(Gnjm) = Qnj " (Gnj) ^ * (R^-njm - RL n j)J • (G n j) 1. (11) 

This equation gives a first-order correction to G n j that is an approximation to G n j m . 

4.2 DETERMINING THE TEMPERATURE DERIVATIVES 
The differential of Eq. (6) yields 

Y £pnj * (dTnm - dTjm)J = dV nm + Y Jcinj ' (dT n m ” dTj m )J (12) 

Where 

^nm = ^Qnm + dDnm + ^ (dac ' Inc) * ^nm * ^n 
+ Y, (dec ' fnc) * An * (Enm _ Bnm) 

- % [K nj * (dT nm - dT jm )] - % [r r ] * (dB nm - dBjm>) (13) 

+ Y (ac * dfnc^ ' ^nm * ^n + ^ ( e c * dfnc) 

• A n • (E nm - B nm ) 

+ Y * ^ nm “ dB nm ) 

where the prefix (d) indicates a differential. The derivative dT nm /dQ nm is obtained 
from Eq. (13) by setting the other differentials to zero. The remaining derivatives are 
obtained by use of the chain rule. For example: 

dTnm/da c = dT nm /dQnm * dQnm/da c (14) 

and 

dQ nm y da c is obtained from Eq. (13). 
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Therefore, all of the sensitivities and all of the derivatives of the temperatures with 
respect to the design variables can be computed with a minimum of work. Some savings 
in effort might be achievable if the derivatives and sensitivities are evaluated only for 
the thermally critical missions; however, there is no clear method for selecting the 
thermally critical missions without first evaluating dT nm /dQ nm . 

The uncertainties in the temperatures are 

<5 ^nm = ^ ^(dT nm /da c • 6a c )2 + ^ (dT nm /de c * 6e c )2^J 1/2 (15) 

where the sums are over the c coatings. The derivatives of the temperature 
uncertainties with respect to the design variables f nc are: 



r dT d 2 T 
run run 

da df da 
L c nc c 


(6 a ) 2 
c 


dT d 2 T 
+ nm run 

de df de 
c nc c 


(6e ) 2 

c 


(16) 


4.3 DETERMINING THE DERIVATIVE OF THE WEIGHT (W) 

The derivatives of W with respect to the fractional area coverages, f nc , determine the 
direction of the step to the next coating pattern (the next design). These derivatives 
are best obtained by the chain rule. The first step is to determine the sensitivities of 
the heater powers to f nc . Suppose we have a set of T nm and 6 T nm for one set of 
heater powers. Suppose the set of T nm and § T nm does not meet the thermal 
requirements. Because we know dT/dQ, we have the approximation: 

^nm/dQnm * (Qnm 1 - Qnm) = ^nm' - Tnm (17) 

where (') indicates the Q or T after the heater power is increased to meet the 
temperature requirements. In this equation, Q and T are known and the temperature 
requirements dictate T'. Therefore, we have a system of linear equations to determine 
Q'. Note that Eq. (17) represents only K equations (one for each temperature-controlled 
node), although these K equations must be solved for each of M missions. Gaussian 
elimination will probably be the most suitable solution method, because there is no 
a priori knowledge of the condition of the matrix, dT/dQ. 
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Solving for Q' would be simple if it were not for the fact that Q' cannot physically be 
negative. If we solve the system of equations for Q' and find that some values are 
negative, we have no method for determining which should be set to zero to minimize 
the weight penalty. This problem is similar to that of the classical regression problem, 
if we are limited to 5 coefficients in a regression involving 10 independent variables, 
which 5 should be chosen? There is no solution to the regression problem, so we can 
expect to find no solution to the Q' problem. Therefore, if Q' is negative we would set 
the most negative Q' to zero. After this is done, we must re-solve the system of 
equations for the remaining Q'. 

The derivative of W with respect to f nc can now be determined. We have from Eq. (2) 

W = max ^ab • Q + bb • £(C n ) + ae • Q + be • £ (C n )J (18) 


so 


dW 

df 


= ab* 


nc 


dQ_ 

df 


aCn 

+ bb *l If - 
nc nc 


+ ae* 


< 10 . 

df 


be*£ 


d C r 


nc 


df 


nc 


The Q in this equation is evaluated for the mission that has the highest total heater 
power. The Cn in this equation is summed over all of the heated nodes (Cn is the 
maximum heater power for node n, considering all of the missions)^ dQ/dfpQ is the 

product of dT/dfnc and dQ/dT from Eq. (17). Similarly, dC n /df nc is determined from 
dQ/dfnc, but for the particular mission that determines each particular heater capacity. 

HA COMPUTING THE NEXT SET OF DESIGN VARIABLES 

The foregoing equations give the objective function and its derivatives. These data 
make hill-climbing techniques, such as the maximum-rate-of-descent method, suitable 
for the optimization. CONMIN (Ref. 5) is a standard computer routine based on hill- 
climbing techniques that is particularly suitable for the thermal optimization problem. 
It is designed especially to handle constraints on the dependent and independent 
variables. 
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The selection of CONMIN followed an analysis of the needs of the thermal optimization 
strategy and a review of CONMIN's capabilities. The requirements were met by 
CONMIN as follows: 


• CONMIN is able to use analytical derivatives as well as finite differences. 

o CONMIN is able to handle non-linear constraints and the particular form of 
the objective function. CONMIN also has a feature that permits the user to 
select the error band on satisfying the constraints, so that the time spent 
satisfying the heat-balance constraint can be controlled. 

e CONMIN has a re-start capability. 

© CONMIN can handle upper and lower bounds on the independent variables. 

• CONMIN has the capability of evaluating all constraints or only those that 
are active. This capability will be especially useful because the tradeoff 
must be made to determine if considering all constraints at every step 
sufficiently reduces the number of erroneous steps. 

© CONMIN has the capability of starting with a nonfeasible solution. It uses a 
penalty function to search for a feasible solution. 

• CONMIN allows for many input constants that permit the user to control the 
parameters of the optimization, such as step size. 

Unlike many optimization problems, the present strategy does not have a problem of 
scaling disparate independent variables, such as pressure and temperature, because all 
of our independent variables are area ratios. 

4.5 ESTIMATED COMPUTATION TIMES 

The practicality of the optimization strategy depends on the cost of performing the 
computations. Two sets of computations are required: the temperatures and the 

temperature derivatives. 

At each step, the program must determine the temperatures for all of the nodes for the 
beginning and end of the mission for M missions. Therefore, the heat-balance equations 
must be solved 2MI times. If there are P models, to incorporate uncertainties in the 
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internal conductances, then the temperatures must be determined 2MPI times. How- 
ever, the way that the powers for K heaters are to be determined requires re-solving 
the heat-balance equations from 1 to K times per set of coatings evaluated. On the 
average, we can expect K/2 re-analyses for the heater sizes. Therefore, the total 
number of solutions to the heat-balance equation will be KMPI for each optimization 
problem. 

Computation of the temperature derivatives requires approximately the same amount 
of work as solving the heat-balance equations, although the equations for the 
derivatives are linear. The derivatives must be determined for each mission, for the 
beginning and end of the mission, and for each model. Therefore, 2MPI computations 
are required. Because we have used the matrix-inversion technique, the sensitivities 
can be determined from the derivatives with little extra work. 

By combining the estimates for the two sets of computations, we find that the problem 
is similar to solving the heat-balance equation MPI • (K + 2) times. For the problem 
size estimated at the beginning of Section 4, the heat-balance equation would be solved 
20 • 4 • 25 * (30 + 2) = 64000 times. This time can be compared to the time now used to 
simulate the transient responses of spacecraft. For example, it is not uncommon to use 
5-minute time steps to analyze the temperature response of spacecraft over 48 hours, 
in search of the quasi-steady-state solution. This computation requires solving the 
heat-balance equations 576 times. If 100 variations are considered, as part of the 
evaluation of the performance for 20 missions, the coating uncertainties, and the 
internal-conductance uncertainties, the optimization time would be comparable to the 
design evaluation time. 

Run times for typical spacecraft for a 48-hour simulation are 1 hour on a VAX1 1/780. 
For the 64000 solutions, approximately 110 hours would be required on the VAX 11/780. 
If a vector processor were used, this time could be reduced by a factor of 10 (Ref. 3). 
Another factor-of-10 time reduction would be realized if a Cyber 175 were used. A 
third factor-of-10 time reduction could be achieved with a Univac 1100/80 with dual 
array processors. On this last machine, the CPU time would be 0.11 hours. As 
indicated in Reference 3, a 190-node problem can be solved 100 times in 8 CPU seconds 
on a CYBER 203. Therefore, 64000 solutions would require 1.42 hours without a vector 
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processor and 0.5 hours with a vector processor. Thus, the run times for the 
optimization will be long but not unreasonable, especially because the output from the 
optimization will reduce the need for some of the studies normally done and because 
the optimization will not be done more than once or twice for an entire project. 
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SECTION 5 


THE OVERALL OPTIMIZATION STRATEGY 

With the mathematics and the computational steps now defined, we can turn our 
attention to the work the designer must do. He must: 

1. Develop thermal models of the vehicle. 

1.1 A baseline model with the nominal values of the conductances 

1.2 Alternative models with off-nominal values 

2. Specify the coatings and coating properties to be used. 

2.1 The nominal values at the beginning and end of the mission 

2.2 The uncertainties at the beginning and end of the mission 

3. Specify the external-surfaces 

3.1 External surface areas 

3.2 Solar and infrared heat fluxes on each surface for each mission 

3.3 The coatings that can be applied to each surface 

4. Specify the temperature-controlled nodes 

4.1 The allowable range in temperature differences 

4.2 The allowable range in the average temperature of groups of nodes 

4.3 The allowable range in differences in average temperatures of groups of 
nodes 

5. Specify the weight penalties for the beginning and the end of of the mission. 

5.1 Allow for decreasing voltage with time due to power-system degrada- 
tion 

5.2 Weight penalty per watt of peak heater power 
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5.3 Weight penalty per watt of average heater power 

6. Specify the optimization parameters for CONMIN 

The maximum step size, the number of steps, etc. 

The coating configuration will be optimized according to the theory described in the 
previous sections and shown schematically in Figure 2 (see also Table II). However, the 
designer uses the optimization procedure only as an aid. He must still exercise his 
judgment: 

• He must select appropriate time periods. In most cases, the time period 
would be one orbit because all sensitive equipment would be well enclosed 
by structure and insulation. Optimization for external equipment, such as 
star trackers and special instruments, might be conducted in a second pass 
using time averages that were close to the time constants of this externally 
mounted equipment. 

• He must select the final heater configuration. The program will assume the 
heaters are applied to every node for which the temperature must be 
controlled. If the heater layout or installation can be simplified, the user 
must perform the simplification. If the heaters cannot be placed on the 
temperature-controlled element (such as a mirror), he must adapt the 
procedure results. For example, he can design a heated plate next to the 
temperature-controlled element. He can then re-run one optimization 
iteration with the heated-plate configuration to test his manual modifica- 
tion. For control of temperature differences, he can imitate what will be 
done in the hardware. For example, if the temperature difference is sensed, 
heat could be applied to the colder node. Alternatively, the designer can 
select the node to which heat would be applied, with the temperature of the 
other node floating. 

• He must optimize the less tangible aspects of the thermal design. For 
example, he must consider the reliability of the design and the appropriate 
design margins, the non-orbital aspects of the missions, design flexibility for 
future changes in requirements or power dissipations, system manufactura- 
bility, system predictability, the costs and schedule requirements for the 
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Figure 2. Flow Chart of Optimization Strategy 
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TABLE II 

POSSIBLE FLOW SCHEMATIC FOR OPTIMIZATION PROCEDURE 


1. Read input data 

• Include all of the missions, alternative models, etc. 

• Include initial values for the coating patterns 

2. Read input data for the CONMIN computations 

• Step size, error bands, etc. 

3. Compute a feasible solution 

• For initial set of design variables, compute the temperatures and heater 
powers for the beginning and end of each mission. Determine the tempera- 
ture uncertainties 

• Compare the temperatures (with uncertainties) to the requirements 

• If solution is not feasible, reduce the solar absorptances and try again 

4. Compute the direction of the next set of design variables: 

• Weight penalty 

• Derivatives of objective function with respect to design variables 

• Derivatives of active constraints with respect to design variables 

• Direction of next step using CONMIN 

5. Compute the new temperatures and heater powers 

• For the new set of design variables, compute the temperatures and heater 

powers for the beginning and end of each mission and determine the 

temperature uncertainties 

e Compare the temperatures (with uncertainties) to the requirements 

• If solution is not feasible, reduce the solar absorptances and try again 

• If solution is feasible, go to Step 4 
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thermal-control equipment and the system development, and the costs of 
testing the thermal-control system. 

e Therefore, the optimization strategy, even when implemented on the 
computer, is a design aid -- not a complete design solution. However, 
because it is a systematic search among alternatives, it is an aid that can 
eliminate many of the intuitive decisions that must be made in the design 
process. 
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SECTION 6 


DEMONSTRATION PROGRAM 


A pilot computer program was developed to validate some of the techniques needed for 
the subsequent optimization: 

I 

• Specifying the coatings and their tolerances 
e Finding a feasible solution 

e Computing the nominal temperatures and uncertainties in temperatures 

• Computing the heater power to exactly meet the temperature requirements 
for all nodes simultaneously 

• Computing the derivatives of the temperatures with respect to the design 
variables. 

Because the purpose of the current work was only to define the optimization strategy, 
the pilot program does not include all of the steps. For example, it does not include an 
algorithm for selecting the next set of design variables. The steps of the logic are 
selected manually, rather than automatically. The program runs interactively. A menu 
appears on the screen and the user selects the next step. The steps include: computing 
the nominal case, computing the effects of uncertainties, selecting the thermostat set 
points, computing the heater-power required and computing the derivatives. 

The pilot program, presented in Appendix A, was written in Microsoft BASIC. It 
probably will require modification for other forms of BASIC, because there is no 
industry standard on this language. 

The pilot program was successfully applied to the three-node sample problem presented 
in Section 3. The details of the application (see Appendix A) confirmed the practicality 
of each step. Although the program tests only selected parts of the computational 
procedure, it does include enough of the steps that we are confident that the entire 
process can be developed and implemented successfully. 
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SECTION 7 


PLAN FOR IMPLEMENTING THE OPTIMIZATION PROCEDURE 


Development of a computer program to implement the entire optimization strategy 
would involve a significant effort. The major steps are defined below. 

1. Design the input and preprocessor to accept a SINDA thermal model (Ref. 4) 
of the vehicle. 

We cite SINDA for use in representing the thermal models because of its 
extensive use in the aerospace industry. Other thermal modeling programs 
could also be used; however, SINDA seems to have significant advantages. 
It is widely used and accepted. It is frequently updated and improved; new 
computational routines can be easily added. And it is applicable to almost 
any thermal problem, including problems with fluid flow (most conveniently 
in a version called SINFLO). On the other hand, if we developed a new 
thermal-analysis subroutine, the user of the optimization process would need 
to learn a new set of input procedures and formats. In addition, he would 
need to have one model for the optimization process and another for the 
detailed SINDA analyses he might conduct later. 

2. Improve the SINDA computational efficiency for multiple-mission analyses. 

SINDA does not have a matrix-inversion solution of the type described in 
Section 4.1. Therefore, such a routine would need to be added. This can be 
accomplished by writing a SINDA subroutine. The capability to accept new 
routines is a standard feature of SINDA. 

3. Compute the derivatives of the temperatures with respect to coating 
properties and coating areas. 
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A tradeoff between numerical and analytical differentiation must be made, 
but we anticipate that the analytical differentiation will be faster. There- 
fore, SIN DA must be modified to compute the derivatives. Some versions of 
SINDA have subroutines for computing sensitivities; these may be suitable. 
A tradeoff is required to determine if the derivatives need to be evaluated 
for all constraints or only for the active constraints. 

4. Modify SINDA to compute the required heater powers. 

A computational subroutine must be added to SINDA for computing the 
required heater powers. Some study is needed to determine which heaters 
should be activated and when all are not needed to meet the constraints. 
The heater power was computed in the pilot program by re-inverting the 
rows of the matrix (Eq. 8) corresponding to the nodes with heaters. A 
similar process can be followed with SINDA, as described in Section 5 ; 
however, the non-linear radiation terms may make this linear technique 
unsuitable. For example, an approximate inverse may be too inaccurate. 
An iteration scheme may be required. 

The computed heater power required to bring the node back into the 
specified range can also be used to compute the gradients needed in 
CONMIN. The matrix of derivatives of the temperatures with respect to 
heater powers can be inverted to give the derivatives of the heater power 
with respect to temperature. Therefore, the temperature derivatives 
obtained in Step 3 can be multiplied by the heater power derivatives to yield 
the heater power derivatives required in the next step. 

5. Modify SINDA to compute the objective function and derivatives of the 
objective function. 

The objective function is computed from Eq. (2). Derivatives are obtained 
from Steps 3 and 4. 
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6. Modify CONMIN to call SINDA for all of the needed information. 

CONMIN, as currently written, would call for the temperatures and the 
derivatives separately. Because many of the computations associated with 
the derivatives are the same as or similar to those associated with the 
temperature computations, CONMIN should probably call SINDA once for all 
of the data: the temperatures (constraints), fraction of the area of each 
external node that is covered by each coating (the design variables), the 
weight penalty (objective function), and the derivatives of the weight 
penalty with respect to the design variables. SINDA will be called twice for 
each mission and for each set of internal conductances so that these values 
can be computed for the beginning and the end of the mission for all possible 
internal conductances. 

The user-selected option in CONMIN to compute all or only the active 
constraints should be utilized. (See the discussion under Step 3.) Note that 
CONMIN permits constraining each design variable independently. Our 
problem also constrains the sum of the fractions of each external area to be 
1.0. It will probably be efficient to analytically eliminate one of the 
independent variables with this constraint. For example, the optimization 
program could be designed to assume that the remainder of the surface is 
coated with a diffuse white surface. 

7. Combine SINDA and CONMIN into a single program. 

The currently separate programs, SINDA and CONMIN, must be combined 
into a single program. It is not clear at present whether SINDA should be a 
subroutine to CONMIN; CONMIN, to SINDA; or both to a new main program. 
We have assumed that a new main program would be used. Much will depend 
on the complexity of the various preprocessors. In any case, the resultant 
combination should probably look like SINDA because the program will be 
used by SINDA users. Extra input data will be required for the CONMIN 
computations. 
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SECTION 8 


PLAN FOR APPLICATION OF PROCEDURE TO SPACE TELESCOPE 


The pilot program gives us considerable confidence that the optimization strategy can 
be implemented in the form of a computer program. Although many of the 
computational steps need development, we show in this section how the computational 
process may be used to optimize the thermal control system of the Space Telescope 
(ST). 

The Space Telescope was chosen for the demonstration because it embodies most of the 
constraints found in a scientific satellite, and its design maturity provides a body of 
data for input to the optimization. The overall plan is to optimize the Space Telescope 
using Perkin Elmer's Space Telescope Systems Thermal Model and a modified set of 
temperature constraints. 

The Space Telescope is a stellar observatory comprised of three major subassemblies: 
the Support Systems Module (SSM), the Optical Telescope Assembly (OTA) and the 
Scientific Instruments (SI) (see Figure 3). There are two major thermal design 
constraints: 


1. maintain all subassemblies within temperature limits that will assure that 
damage does not occur due to high or low temperatures, and 

2. maintain the change in temperatures between components, within a sub- 
assembly, within limits that assure optical alignment during any twenty-four 
hour period. 

Maximum/minimum temperature limits must be maintained over all phases of the ST 
mission (launch, deployment, scientific operation, etc.) while optical alignment must 
only be maintained during the scientific operations which occur when the ST is in the 
orbit shown in Figure 4. During scientific operations the solar line of sight can fall 
within the 135 degree envelope shown in Figure 4. As the vehicle moves relative to the 
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Figure k. Space Telescope Orbit Parameters 



sun, the change in solar viewing angle will cause the Space Telescope to heat or cool by 
as much as 20OF. Since the major driver in any design is the scientific portion of the 
mission, the temperature level design should only be optimized for the operational 
phase of the mission. The resulting design should then be compared to the launch and 
deployment requirements so that the specific design changes required to meet the 
nonoperational constraints can be developed. 

The SSM provides the overall enclosure, power, power control, guidance and position 
stability. Its thermal control requirements and nodal breakdown are presented in 
Figures 5 and 6. Because it does not provide a structure for any of the optical 
components, the allowed temperature ranges and levels for the SSM are significantly 
greater than those for either the OTA or the Sis. The OTA provides both the structure 
to hold the optical components (i.e., primary mirror, etc.) and the Sis. Because of the 
optical requirements the allowed temperature variations are small. These requirements 
also tend to limit the choice of materials to those which have a low thermal coefficient 
of expansion (Graphite Epoxy) and sometimes require the use of resistance heaters to 
maintain alignment. These all tend to minimize a component's allowed operational 
temperature range. The SI provides the structure that holds the experiment and its 
required optics and electronics. While not as massive as the OTA its structural and 
resulting thermal requirements are similar to the OTA'S. Both the thermal require- 
ments and a nodal breakdown for the the OTA and a typical SI are given in Figures 7 
through 12. 

A significant aspect of the optimization problem is the uncertainty in the specified 
coating properties. This is included by utilizing the thermal properties and the 
uncertainties presented in Table III. 

The optimization strategy is to minimize the spacecraft's weight (Eq. 2). For cold bias 
systems that are maintained by resistance heaters, this requires knowledge of the 
relationship of weight to power (see Table IV). Three relationships derived from ST 
data are used to determine the weighting factors ab, ae, bb, and be: 

• AVERAGE POWER - This is the penalty for the average power required to 
maintain the spacecraft. The size of the solar array, battery, power 
distribution and control etc., weights are accounted for in this term. 
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Figure 5. Space Telescope Support Systems Module (SSM) Thermal Node Definition 
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Figure 7. Space Telescope, Optical Telescope Assembly - Forward Section Thermal Node Description 
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Figure 8. Space Telescope, Optical Telescope Assembly - Central and Main Baffle Thermal Node Description 
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Figure 9. Space Telescope, Optical Telescope Assembly - Metering Truss and Secondary Mirror Thermal Node 
Description 
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Capacitance 

BTU/°F 

Surface 

Area 

FT 2 

Temperature 

Requirements 

°F 

Surface 

Inner 

Finish 

Outer 

Power 

Dissipation 

Typical Axial^^ 
SI Module 
Optical Bench 

3011 

74 

111 

X±1 

None 

MLI2 

150/5S< 2 > 

Radiator 

3012 

18 

44 

None 

Black Paint 

Black Paint 

None 

Interior 

3013 

16 

44 

None 

Black Paint 

MLI2 

None 

+V[ Ends 

3014 

6 

11 

None 

Black Paint 

MU2 

None 

"B" Interface 

3015 

P 

0 

i None 

None 

None 

None 


NOTE: (1) One of eight axial and radial modules. 


(2) The optical bench power dissipation will range from 58 to 150 Watts. 


Figure 10. Space Telescope, Optical Telescope Assembly - Scientific Instruments Thermal Node Description 
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(NO THERMAL CONTROL) 
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Thermal 

Surface 
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Node 

Capacitance 

Area 

Requirements 

Surface 

Finish 

Power 

Component 

Number 

BTU/°F 

FT 2 

°F 

Inner 

Outer 
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+V3 Truss 

V3 Barrier 

-V3 Truss 
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Shelf 

Diamond Truss 

Diamond Truss 
Interface 


Figure 12. Space Telescope, Optical Telescope Assembly - Aft Focal Plane Thermal Node Description 















TABLE HI 

PROPERTIES OF AVAILABLE THERMAL COATINGS 





Alpha 


Epsilon 

No. 

Description 

Mean 

Variation 

Mean 

Variation 

1 

FOSR 

0.1 

+0.05 

0.85 

+0.05 

2 

Aluminum 

0.2 

+0.05 

0.10 

+0.01 

3 

Max Orb 

0.95 

+0.05 

0.1 

+0.01 

4 

Black Paint 

0.95 

+0.02 

0.95 

+0.02 

5 

Gold 

0.25 

+0.03 

0.03 

+0.01 

6 

Multi-Layer 
Insulation (high a) 

0.95 

+0.05 

0.01 

+0.0 

-0.005 

7 

Multi-Layer 
Insulation (low a) 

0.1 

+0.01 

0.01 

+0.000 

-0.005 

8 

Optical Surface 

0.10 

+0.02 

0.02 

+0.01 

9 

Glass 

0.1 

+0.05 

0.9 

+0.05 
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TABLE IV 

POWER PENALTIES 


Parameter 

Penalties 


Peak Power 

0.15 

lbm/watt 

Average Power 

0.7 

lbm/watt 

Life Cycled) 

0.56 

watts end of mission 



watts start of mission 


Based on a 32 to 2k supply voltage variation (start to end) and 
resistance heaters 


( 1 ) 


• PEAK POWER - This is the penalty for the peak power required to maintain 
the spacecraft. The number of heaters and their duty cycle has a significant 
bearing on the weight of a thermal control system and this is accounted for 
in this term. 

• LIFE CYCLE - As the solar arrays age their output and the corresponding 
supply voltage decreases. This results in a change in both the average and 
the peak power availability with time. The weight penalty for power vs 
spacecraft life is accounted for in this term. 

The factors given above are combined as shown below to determine the power weighting 
factors: 


ab = 

average power 

ae = 

average power/life cycle 

bb = 

peak power 

be = 

peak power/life cycle 


After the model and weighting factors are defined, the model is put in the form to be 
optimized as given in Table V. 

The ST will be optimized for two missions, each with the solar vector at opposite ends 
of the allowable envelope shown in Figure 4. During this process the coating on the 
SSM will be optimized to meet the temperature requirements listed in Figures 5 through 
10. If the absolute level is not specified the upper and lower bounds for these 
components will be set at 70°F plus or minus the prescribed temperature tolerance. 
For example, the limits for node five in Figure 7 are: 

L 5 = 70 -0.6 = 69.4°F 

U 5 = 70 +0.6 = 70.6°F 

When the analysis for the average temperature of 70°F is complete, other average 
component temperatures will be tried to assess the impact of temperature level on the 
optimal design. During this process we will not vary the coatings on nodes other than 
the exterior surface since these are generally set by requirements other than thermal. 
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TABLE V 

OPTIMIZATION MODEL 


• Model size 

e Nodes for coating optimization 

(f nc varied) 

« Nodes with temperature 
requirement, constraints 


• Nodes with weighted - average 
temperature constraints 


• Nodes with temperature 
difference constraints 


• Nodes with a weighted average 
temperature difference 
constraints 


84 Nodes 

7 SSM Nodes (Figure 5) 

29 a OTA and SI Nodes 

none 

2 

difference between the metering truss 
(Mid) and metering truss (Fwd) 

none 


a The twenty-nine nodes requiring temperature level control are shown in 
Figures 7 through 12. The optical benches of the 3 axial SI modules and 4 
radial modules (not shown) also require thermal control. 
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However, the program will determine the optimal placement of heaters, and their 
location will not be fixed by their locations in the ST design. 

REMARKS ON ST THERMAL DESIGN PROBLEM 

During its design, the ST spacecraft was partitioned into three major subsystems. Each 
subsystem's contractor tried to optimize that subsystem within that subsystem's 
constraints. For example, the SSM designer had a wide latitude in temperature 
excursion during an observational period and over the entire life of the mission. This 
led to a passive SSM shell design. On the other hand, faced with a varying environment, 
both the OTA designer and the SI designer were led to active systems to maintain both 
level and change over any period. Additionally, because of the problems of maintaining 
the figure of a mirror at temperatures other than the temperature at which it was 
figured, a significant portion of the OTA is maintained at a single temperature of 70°F. 
Optimization could change parts of this. Since the entire ST would be optimized as an 
entity, it is possible that maintaining the SSM at a fixed temperature (plus control 
tolerance) may be a lower weight solution than the current design. Additionally it 
would be of interest to modify the OTA criteria from one that requires a 70°F structure 
to one that only maintains the mirrors at that temperature and finds the optimum 
temperature for the remaining structure. This may result in a telescope with a 70°F 
Primary Mirror and a different temperature Focal Plane Assembly. 

The accepted thermal control procedure is to find the maximum temperature of a 
component and then maintain it using heaters and a single set point thermistor 
controller. This assures that all temperature and gradient constraints are met for all 
phases of the mission. However a more power-efficient method of control may be to 
maintain the temperature difference between components by using the temperature of 
one component as the set point temperature of the other. This type of control is 
possible with presently available micro-computers and may result in a significant power 
(thus weight) savings for the optimized spacecraft. 
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SECTION 9 


CONCLUSION 


The purpose of the work reported here was to develop a strategy by which thermal 
designers for spacecraft could devise an optimal thermal control system to maintain the 
required temperatures, temperature differences, changes in temperatures and changes 
in temperature differences for specified equipment and elements of the spacecraft's 
structure. Thermal control would be maintained by the optical coating pattern chosen 
for the external surfaces and the heaters chosen to supplement the coatings. A 
strategy appropriate to computer-aided design was anticipated. 

We focused on the development of a practical tool for the designer: an optimization 
program compatible with his present analytical tools and useful in the conceptual- 
design phase. We initially examined the general problem of spacecraft temperature 
control, including all phases of the mission and all phases of the design process. To aide 
in conceptualizing the problem, we set up a simple three node spacecraft to optimize. 
This provided us with insight into the problem without the confusion of a complex 
model. From this exercise we learned that two major problems needed to be overcome: 
(1) how to optimize both power and weight, and (2) how to include coating uncertainties. 
These were solved by expressing both power and weight as cost functions to be 
minimized and including both variable properties and uncertainties in the properties of 
the coating. 

We then reviewed the methods available for computation. The optimization requires a 
considerable amount of computer time to solve a complex problem. However, with 
estimated run times of 0.5 to 10 hours, the solution times are not considered a 
significant driver in the problem. This is especially true considering the potential for 
savings in the development of a space system when optimization is used. 

The next task was to develop the overall optimization strategy. The strategy developed 
does not eliminate the need for a thermal designer, but is a design tool that eliminates 
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many of the intuitive decisions that must be made. While the work to this point 
demonstrated the feasibility of the proposed optimization procedure, a pilot program 
was developed to validate some of the techniques needed for the optimization. The 
pilot program was successfully applied to the three-node sample problem presented in 
Section 3. The details of the application (see Appendix ) confirmed the practicality of 
each step. Although the program tests only selected parts of the computational 
procedure, it does include enough of the steps so that we are confident that the entire 
algorithm can be developed. We performed one additional test on the optimization 
process by conceptually applying it to the Space Telescope. The work remaining is 
detailed in Section 7. 

In developing the optimization strategy, we have included all of the elements normally 
considered by the designer: vehicle orientations, equipment power dissipations, 

coating-property uncertainties, uncertainties in the internal thermal conductances, and 
coating-property changes. The method uses the same thermal model data the designer 
develops for his analyses; there is no extra work. He can use his intuition to speed the 
optimization process and to simplify the resultant heater arrangement, but he can also 
use the process as a "black box." We have demonstrated how the optimization process 
can be cast in computer form; the results will be in the form of a specification of the 
coating pattern to be used on each external surface, the heater capacities required for 
each node, the average power required for each mission, and the maximum and 
minimum temperatures expected for each node. The method is limited, however, to 
quasi-steady-state temperatures; therefore, the designer must select time periods over 
which steady-state temperatures are representative of the mission temperatures. 
Because most sensitive equipment is located well inside the vehicle, his selection is 
usually not difficult: the orbital period is adequate. 

Three of the four steps to having an operational optimization program have been 
accomplished: 

o The governing equations have been developed 
e The theory has been developed for solving the equations 
o The methods have been demonstrated. 
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The last but largest step - developing the production program, including making the 
final tradeoffs and determining how to integrate the theory with existing analysis and 
optimization programs - still must be accomplished. 
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SECTION 10 


NOMENCLATURE 

A surface area 

B black body radiation factors 

C heater capacity 

D electronic power 

E earth flux 

G weight average temperature 

^nj K n j + RL n j 

K thermal conductance 

K n j conductance between nodes n and j 

L lower temperature limit 

M mission 

N total number of nodes 

P model version 

Q the heater energy required for the most severe environment; maximum 

total power demand 

Qnm orbital average power demand for node (n) in mission (m) 

R radiation interchange factor 

RLnjm linearized radiation term (Eq. 8) evaluated for temperatures in mission 
m 

R^-nj Linearized radiation term (Eq. 8) evaluated for maxiroum temperatures 
S solar flux 
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NOMENCLATURE (Continued) 


Tj m temperature of node j in mission m 

T nm temperature of node n in mission m 

U upper temperature 

V nm residual in heat balance equation 

W weight 

a solar absorptance 

ab weight per unit of power, beginning of life 

ae weight per unit of power, end of life 

bb weight per unit of power (capacity) beginning of life 

be weight per unit of power (capacity) end of life 

c coating 

ca changes in absorptance, beginning to end of mission 

ce changes in emissivity, beginning to end of mission 

c n heater capacity of node (n) 

da range of uncertainty in solar absorptance, beginning to end of mission 

de range of uncertainty in emissivity, beginning to end of mission 

6a c uncertainty in solar absorptance 

Se c uncertainty in emittance 

6T nm uncertainty in temperature 

e emissivity; infrared emittance 

f nc fractions of the area of node (n) that are covered with coating (c) 
j node j 

m mission number 

n node n 
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APPENDIX A 


PILOT PROGRAM 


In developing the optimization strategy presented in Sections 3 and 4, we were careful 
to test the practicality of each step. The testing was done with a pilot program that 
used selected parts of the technique. In particular, we selected those parts that 
demonstrated computing the nominal case, computing the effects of uncertainties, 
selecting the thermostat set points, computing the heater-power required, and 
computing the derivatives. 

The pilot program was tested with the three-node satellite used as an example in 
Section 3 (see Fig. A-l). Two of the nodes are external, one facing toward the earth 
and the other facing away from the earth. The third node is internal. Each node is 
connected to the other two. To minimize computation time, all of the radiation terms 
were linearized. 



Figure A-l. Three-Node Problem 
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The input data for the three-node, vehicle are listed in Table A-l. In response to 
prompts from the program, the user specifies the number of nodes, number of missions, 
and the coating properties. The coating properties were assumed to be known to within 
plus or minus 0.03. The user also specifies the areas of the external nodes and the 
starting values of the f nc 's, the fractions of external node areas covered by each of the 
coatings. The input calls for an initial guess at the temperature, in anticipation of an 
iterative solution technique. Because a matrix-inversion technique is used, this initial 
temperature does not enter the solution for the pilot program. The internode 
conductances are assumed to be symmetric. The missions are described by the 
electrical heat input to each node from the electronic equipment, the solar flux 
impinging on each node, and the infrared flux emitted by the earth that impinges on 
each node. Only two missions were included in the sample case. All of the values used 
in the sample case were the default values; therefore, there are no entries in response 
to the prompt asking if there is a change to be made. The computed conductance 
matrix is listed at the bottom of the input. 

The output data for the three-node vehicle are listed in Table A-2. As a check, the 
residual of the conductance matrix after inversion and the inverted matrix are listed. 
The residual is seen to be the identity matrix, so the error in the inversion process was 
insignificant. The pilot program calculates the nominal temperatures with no heater 
input immediately. The temperatures are listed after the inverted matrix. For 
example, the temperature of the internal node, Node 3, is 64.7°F for Mission 1. The 
chart of temperature data includes the area ratios for the various coatings and the 
effective solar absorptance and an infrared emittance (0.337 and 0.421 , respectively, 
for Node 1). Because there are no heaters for this nominal case, there is no weight 
penalty, as indicated in the output below the temperature chart. 

In response to a list of options displayed on the CRT screen, but not shown in the output 
data, we chose to compute next the effects of the coating uncertainties (called 
tolerances in the output data). The temperatures with the uncertainties included are 
tabulated under the heading, "IF COATING TOLERANCES ARE INCLUDED." The MIN 
temperature, for example, is the nominal temperature less the three-sigma temperature 
uncertainty. For Node 3, the temperatures will be between 38.3°F and 9i.l°F in 
Mission 1 and between 38.6°F and 87.90F in Mission 2, if no heaters are used and this 
coating pattern is chosen. 
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TABLE A-l 

INPUT DATA FOR THE DEMONSTRATION CASE 


NUMBER OF NODES 3 CHANGE TO? 

YOU CAN HAVE UP TO 51.1404 MISSIONS 
NUMBER OF MISSIONS: 2 CHANGE TO? 


THE FOLLOWING COATINGS ARE AVAILABLE: 


NO. 

DESCRIPTION 

ALPHA 

. EPSILON 

1 

FOSR 

. 1 

.85 

2 

ALUMINUM 

.2 

.1 

-» 

MAXORB 

.95 

. 1 

4 

BLACK PAINT 

.95 

.95 


NODE PROPERTIES: 

FOR NODE NO. 1 : 

EXTERNAL AREA* A: 300 CHANGE TO? 

INITIAL TEMPERATURE* • T IN DEG F: 70 CH M »^i TO? 

AREA FRACTION WITH AL= .1 *EP= .85: .36 CHANGE TO? 

AREA FRACTION WITH AL= .2 * EP= .1 : .41 CHANGE TO? 

AREA FRACTION WITH AL= .95 *EP= .1 : .17 CHANGE TO? 

AREA FRACTION WITH AL= .95 * EP= .95 s .06 CHANGE TO? 

FOR NODE NO. 2 : 

EXTERNAL AREA* A: 300 CHANGE TO? 

INITIAL TEMPERATURE* T IN DEG F: 70 CHANGE TO? 

AREA FRACTION WITH AL= .1 * EP= .85 : .36 CHANGE TO? 

AREA FRACTION WITH AL= .2 *EP= .1 : .32 CHANGE TO? 

AREA FRACTION WITH AL= .95 *EP= .1 : .17 CHANGE TO? 

AREA FRACTION WITH AL= .95 *EP= . 95 : .15 CHANGE TO? 

FOR NODE NO. 3 : 

EXTERNAL AREA* A: 0 CHANGE TO? 

INITIAL TEMPERATURE. T IN DEG F: 70 CHANGE TO? 

CONDUCTION MATRIX: 

C( 1 * 2 ) 

C( 1 * 3 ) 

C< 2 * 3 ) 

MISSION CHARACTERISTICS: 

MISSION NODE ELEC INPUT 

BTU/HR 

11 0 

12 0 

1 3 0 

MISSION NODE ELEC INPUT 

BTU/HR 

2 10 

2 2 0 

2 3 0 

COMBINED CONDUCTANCE MATRIX: 

-281.58 125.00 125.00 

125.00 -2B7.31 125.00 

125.00 125.00 -250.00 


CHANGE TO? 
CHANGE TO? 
CHANGE TO? 


EARTH FLUX 
BTU/HR-SF 

7.33333 CHANGE (Y/N)? 
41 CHANGE < Y/N ) ? 

0 CHANGE (Y/N)? 

EARTH FLUX 
BTU/HR-SF 

23.3333 CHANGE (Y/N)? 
23.3333 CHANGE (Y/N)? 
0 CHANGE < Y/N ) ? 
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TABLE A-2 

OUTPUT DATA FOR THE DEMONSTRATION CASE 


5I!P UAL IN THE THAT WAS INVERTED: 

1.000000 0.000000 0.000000 

0.000000 1.000000 0.000000 

0.000000 0.000000 1.000000 


INVERSE MATRIX: 

-0.015950 -0.013303 -0.014627 
-0.013303 -0.015543 -0.014423 
-0.014627 -0.014423 -0.018525 


MISSION: 1 

NODE TEMPERATURE 
MIN (F) MAX 

1 92.4 92.4 

2 36.9 36.9 

3 64.7 64.7 


HEATER PWR 

HEATER CAP 

ALPH 

EP 

AREA 

RATIO 

BY COATING NO 

BTU/HR 

0 

0 

0 

BTU/HR 

0 

• 337 

• 421 

1 

0-360 

2 

0.410 

3 4 

0. 170 0.060 

0 

0 

.404 

• 498 

0.360 

0.320 

0. 170 0. 150 


MISSION: 2 

NODE TEMPERATURE 
MIN (F) MAX 

1 63.0 63.0 

2 63.6 63.6 

3 63.3 63.3 


HEATER PWR 
BTU/HR 


0 

0 

0 


HEATER CAP 
BTU/HR 


0 

0 

0 


ALPH EP 

.337 .421 
.404 .498 


AREA RATIO BY COATING NO. 
12 3 4 

0.360 0.410 0.170 0.060 
0.360 0.320 0.170 0.150 


PREVIOUS WEIGHT FUNCTION WAS 0 
THE NEW WEIGHT FUNCTION IS 0 


BTU/HR EGUIV AVERAGE HEATER POWER 


MISSION: 1 

NODE TEMPERATURE 
MIN (F) MAX 

1 64.4 120.5 

2 12.2 61.7 

3 38.3 91.1 


IF COATING TOLERANCES ARE INCLUDED 


HEATER PWR 

HEATER CAP 

ALPH 

BTU/HR 

BTU/HR 


0 

0 

.337 

0 

0 

.404 

0 

0 



EP AREA RATIO BY COATING NO. 

1 2 3 4 

•421 0.360 0.410 0.170 0.060 

.498 0.360 0.320 0.170 0.150 


MISSION: 2 

NODE TEMPERATURE 
MIN (F) MAX 

1 38.1 87.9 

2 39.2 88.0 

3 38.6 87.9 


HEATER PWR 
BTU/HR 

0 

0 

0 


HEATER CAP 
BTU/HR 

0 

0 

0 


ALPH EP 

.337 .421 
.404 .498 


AREA RATIO BY COATING NO. 
12 3 4 

0.360 0.410 0.170 0.060 
0.360 0.320 0.170 0.150 


PREVIOUS WEIGHT FUNCTION WAS 0 
THE NEW WEIGHT FUNCTION IS 0 


BTU/HR EQUIV AVERAGE HEATER POWER 


ENTER 

the thermostat 

SET 

POINTS: 

NODE 

1 

SET POINT 
-999 

(F) 

CHANGE 9 ™^ N ° THERM0STAT 0R HEATER) 

2 

-999 


CHANGE TO? 

3 

-999 


CHANGE TO? 70 
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In running the demonstration case, we next wanted to determine the heater power 
required to hold the internal temperature (Node 3) to 70°F. The thermostat set points 
were entered and the program calculated the required heater power (including the 
effects of the uncertainties) and the resultant nominal temperatures. Notice that the 
heater power for Mission 1 is 1710.42 Btu/hr. If this power were dissipated and the 
coating properties were their nominal values, Node 3 would be at 96.40F for Mission 1. 
Of course, the thermostatic control would take effect, so this much power would not be 
dissipated. From the previous data, we saw that with no heat input, Node 3 would be at 
64.7°F with the nominal coating properties. Therefore, if the coating properties were 
the nominal values, a slight amount of heat would be required to maintain 70°F. 
However, we also saw that, with no heater power and the upper limit of the coating 
uncertainties, Node 3 could be as high as 91. IF. If this were unacceptable, a new set of 
coatings would be selected. 

If we return to the system performance with the heaters on, but we add the coating 
uncertainties (tolerances), the last set of temperature output shows that 69.5°F can be 
maintained at Node 3 for Mission 1 and 69.4°F, for Mission 2. These results are not 
exactly 70°F because the pilot program does not include the method presented in 
Section 4 by which the sensitivity of the temperature uncertainty to the coating 
properties can be computed. 

Finally, in running the pilot program, we selected the option of having the derivatives 
calculated, as shown in the last set of output data. These are the derivatives of the 
temperatures with respect to the f nc 's, the fractions of nodes 1 and 2 that are covered 
by each coating. These derivatives could be combined with the derivatives of the heater 
powers with respect to temperature (dq/dt) and the derivatives of the weight penalty 

with respect to the heater powers (dw/dq) to estimate the next best set of fnc' s * 
However this logic was not included in the pilot program. 

Again, the purpose of the pilot program was not to demonstrate the complete logic but 
to test the ideas. It kept the development of the optimization strategy firmly 
implanted in reality. The success of the pilot program is indicative of the success we 
can expect in the total implementation of the method. 

The complete listing of the pilot program, in Microsoft BASIC, is presented in 
Table A-3. 
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TABLE A-3 


THERMAL OPTIMIZATION PILOT PROGRAM 


00010 ’PROGRAM EXAMPLE/OPT ON PERKINEL AND MINDISK2 

00020 ' JUNE 2, 1983 8:52 PM 

00030 'THIS PROGRAM TAKES A SMALL SATELLITE THERMAL MODEL AND OPTIMIZES THE 
00040 'ALPHA, EPSILON AND HEATER CAPACITY FOR THE NODES SELECTED. 

00050 'TO INVERT A 10-NODE MATRIX REQUIRES 27 SECONDS 
00060 DEFINT I-N 

00070 L$="OFF" 

00080 INPUT "DO YOU WANT TO USE AN EXISTING DATA FILE (Y/N)";SI$ 

00090 IF SI$ "Y" THEN SI$="N":GOTO 130ELSE INPUT "WHAT IS THE FILE NAME";SI$ 
00100 OPEN "I",1,SI$ 

00110 INPUT#1,N,M,NC 

00120 NS=N:MS=M 

00130 INPUT "WILL YOU WANT TO SAVE THE INPUT DATA ON DISK (Y/N)";S$ 

00140 IF S$ "Y" THEN S$="N":GOTO 150ELSE INPUT "UNDER WHAT FILE NAME";S$ 
00150 C$="N": INPUT "DO YOU WANT TO TURN THE LINE PRINTER ON (Y/N)";C$:PRINT 

00160 IF C$ "Y" THEN 190 

00170 INPUT "WHAT TITLE DO YOU WANT";C$ 

00180 PRINT CHR$ ( 1 4 ) ; TAB (40-LEN (C$)/2);C$: L$="ON" : PRINT : PRINT 
00190 IF SI$="N" THEN N=3:M=2:NC=4 'DEFAULT VALUES 

00200 N1=45:N2=55:CH$="CHANGE TO" 

00210 PRINT "NUMBER OF.NODES";TAB(Nl);N;TAB(N2);CH$; sINPUT N 

00220 IF N= NS THEN 240 

00230 PRINT"CAN'T HAVE NEW CASE WITH FEWER NODES THAN STORED FILE:"NS:GOTO 210 
00240 DIM C(N,N) ,A(N) ,AL(NC) ,EP(NC) ,NC$(NC) ,H(N) ,CI(N,N) ,SP(N) ,HC(N) 

00250 DIM CT(N,N),CN(N,N),AR(N,NC) 

00255 1=2975 

00260 PRINT TAB(10);"YOU CAN HAVE UP TO"; (I-4*N*N-8*N)/(4*N*N+7*N); "MISSIONS" 

00270 PRINT "NUMBER OF MISSIONS :";TAB(N1) ;M;TAB(N2);CH$; sINPUT M 

00280 IF M= MS THEN 300 

00290 PRINT"CAN'T HAVE NEW CASE WITH FEWER MISS'S THAN STORED FILE:"MS:GOTO 2 

00300 DIM T(N,M),Q(N,M),S(N,M),E(N,M),QH(N,M),NH(N,M),TU(N,M),DE(N,M,N,NC) 

00305 'THE ARRAYS TAKE 8*N+(4*M+4)*N*N+7*N*M STORAGE LOCATIONS 

00310 IF SI$="N" THEN 400 

00320 'READ THE DATA FROM THE DISK 

00330 FOR 1=1 TO NC:INPUT#1,NC$(I),AL(I),EP(I):NEXT 

00340 FOR 1=1 TO NS-l:FOR J=I+1 TO NS:INPUT#1,C(I,J):NEXT:NEXT 

00350 FOR 1=1 TO NS:INPUT#1,A(I),T(I,1):IF A(I)=0 THEN 370 

00360 FOR J=1 TO NC:INPUT*1,AR(I,J):NEXT 

00370 NEXT 

00380 FOR 1=1 TO MS:FOR J=1 TO NS:INPUT#1 ,Q(J,I) ,S(J,I) ,E(J,I):NEXT:NEXT 

00390 FOR 1=1 TO NS:INPUT#l,SP(I):NEXT:CLOSE:GOTO 550 

00400 IF N 3 OR M 2 THEN 550 

00410 FOR J=1 TO NC:READ NC$(J),AL(J),EP(J):NEXT 

00420 DATA FOSR.O. 1 ,0.85, ALUM, 0.2, 0.1 .MAXORB, 0.95, 0.1 .BLACK PAINT, 0.95, 0.95 
00430 C(1,2)=125:C(1,3)=125:C(2,3)=125 'DEFAULT VALUES 

00440 T( 1 , 1 )=70:T(2, 1 )=70:T(3, 1 )=70 

00450 AR( 1 , 1)“.36:AR(1 ,2)=.41:AR( 1 ,3)=. 17:AR(1 ,4)=.06 

00451 AR(2,1)=.36:AR(2,2)=.32:AR(2,3)=.17:AR(2,4)=.15 

00460 A(1)=300:A(2)=300 

00470 S(1,1)=80000/300:S (2, 1)=5000/300:E(1,1)=2200/300:E( 2,1 >=12300/300 
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00480 S(1,2)=40000/300:S(2,2)=40000/300:E(1,2)=7000/300:E(2,2)=7000/300 

00490 FOR 1-1 TO N:SP(I)=-999+460:HC(I)=0:FOR J=1 TO M:QH(I , J)=0:NEXT:NEXT 
00500 *1. ENTER THE PROPERTIES OF THE EXTERNAL SURFACES 

00510 PRINT: PRINT "THE FOLLOWING COATINGS ARE AVAILABLE:" 

00520 PRINT " NO. ";TAB ( 5) ; "DESCRIPTION" ; TAB (20) ; "ALPHA" ; TAB (30) ; "EPSILON" 
00530 FOR 1=1 TO NC:PRINT TAB(2) ; I;TAB(7) ;NC$(I) J TAB (20) ;AL(I) ;TAB(30) ;EP(I):NE 
00540 PRINT 

00550 PRINT: PRINT "NODE PROPERTIES:" 

00560 FOR 1=1 TO N 

00570 PRINT "FOR NODE N0.";I;":" 

00580 PRINT TAB( 5); "EXTERNAL AREA, A: ";TAB (N1 ) ; A(I) ;TAB(N2) ; CH$ ; : INPUT A(I) 
00590 PRINT TAB(5);"INITIAL TEMPERATURE, T DEG F:";TAB(N1);T(I,1);TAB(N2);CH 
00600 INPUT T(I , 1 ) :IF A(I)=0 THEN 630ELSE FOR J=1 TO NC 

00610 PRINT TAB ( 5 ) J "AREA FRACTION WITH AL=" ; AL(J) ; ",EP=" ;EP( J) ;":";TAB(N1 ) ; 
00620 PRINT AR(I,J);TAB(N2);CH$; : INPUT AR(I,J):NEXT 
00630 NEXT 

00640 FOR 1=1 TO N:F0R J=1 TO NC:IF A(I)*AR(I,J) 0 THEN 680ELSE NEXT J:NEXT I 
00650 PRINT "FOR THERE TO BE A SOLUTION, AT LEAST ONE"; 

00660 PRINT " NODE MUST HAVE AN EXTERNAL AREA.": GOTO 560 

00670 '2. ENTER THE CONDUCTION MATRIX 

00680 PRINT: PRINT "CONDUCTION MATRIX:" 

00690 FOR 1=1 TO N-1:F0R J=I+1 TO N 

00700 PRINT TAB ( 20 ) ; "C ( " ; I ; " , " ; J ; " ) = ";TAB(N1);C(I,J);TAB(N2);CH$;:INPUT C(I,J 
00710 C(J,I)=C(I,J) 

00720 NEXT: NEXT 

00730 '3. ENTER THE INTERNAL HEATING RATES AND ENVIRONMENTAL FLUXES 

00740 PRINT: PRINT "MISSION CHARACTERISTICS :":F0R 1=1 TO M 

00750 PRINT "MISS ION"; TAB ( 10) ; "NODE" ; TAB (20) ; "ELEC INPUT" ;TAB( 35); "SOLAR FLUX"; 
00760 PRINT TAB (50); "EARTH FLUX" 

00770 PRINT TAB ( 20) ; "BTU/HR" ;TAB (35) ; "BTU/HR-SF" ; TAB ( 50) ; "BTU/HR-SF" 

00780 FOR J=1 TO N 

00790 PRINT I;TAB(10) ; J;TAB(20) ;Q(J,I) ;TAB(35) ;S(J,I) ;TAB(50) ;E(J ,1) ;TAB(60) ; 
00800 C$="N": INPUT "CHANGE (Y/N)";C$:IF C$="N" THEN 820 

00810 PRINT TAB (20);: INPUT Q(J,I) ,S(J,I) ,E(J,I) 

00820 NEXT: NEXT 

00830 '4. SET UP THE CONDUCTANCE MATRIX 

00840 FOR 1=1 TO N:C(I,I)=0:FOR J=1 TO N:IF J=I THEN 860 

00850 C(I,I)=C(I,I)-C(I,J) 

00860 NEXT J 

00870 H(I)=.001713*((T(I,l)+459.67)/100) 3 

00880 H(I)=.25 'TEMPORARY FOR CHECKOUT 

00890 NEXT I 

00900 ’5. FOR THE GIVEN AL,EP, FINISH THE CONDUCTANCE MATRIX 

00910 FOR 1=1 TO N:FOR J=1 TO NC:C(I,I)=C(I,I)-H(I)*A(I)*AR(I,J)*EP(J):NEXT:NEX 

00920 FOR 1=1 TO N:FOR J=1 TO N:CT(I,J)=C(I,J):NEXT:NEXT 

00930 PRINT: PRINT "COMBINED CONDUCTANCE MATRIX:" 

00940 FOR 1=1 TO N:FOR J=1 TO N: PRINT USING "#######.##" ;C (I , J) ; : NEXT: PRINT :NEX 
00950 '6. INVERT THE CONDUCTANCE MATRIX 

00960 FOR 1=1 TO N-l:FOR J=I+1 TO N:CI(I,J)=0:CI(J,I)=0:NEXT:CI(I,I)=1:NEXT 
00970 CI(N,N)=1 

00980 '6.1. FIND THE LARGEST DIAGONAL 

00990 CX=0:FOR 1=1 TO N:IF CT(I,I) CX THEN CX=CT(I ,1) :K( 1 )=I 
01000 NEXT 

01010 '6.2. PERFORM THE INVERSION USING THE LARGEST DIAGONAL FIRST 
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01020 

01030 

01040 

01050 

01060 

01070 

01080 

01090 

01100 

OHIO 

01120 

01130 

01140 

01150 

01160 

01170 

01180 

01190 

01200 

01210 

01220 

01230 

01240 

01250 

01260 

01270 

01280 

01290 

01300 

01310 

01320 

01330 

01340 

01350 

01360 

01370 

01380 

01390 

01400 

01410 

01420 

01430 

01440 

01450 

01460 

01470 

01480 

01490 

01500 

01510 

01520 

01530 

01540 

01550 


'DEFER SEARCH FOR LARGEST DIAGONAL TO A LATER VERSION OF THIS PROGRAM 
FOR K=1 TO N-l :FOR I=K+1 TO N 

R=CT (l,K)/CT(K,K): FOR J=1 TO N:CT(l,J)=CT(I,J)-R*CT(K , J) 

CI(I, J)=CI(I,J)-R*CI(K, J) 

NEXT: NEXT: NEXT 

FOR K=N TO 2 STEP -l:FOR I=K-1 TO 1 STEP -1:R=CT(I,K)/CT(K,K) 

FOR J=N TO 1 STEP -l:CT(l, J)=CT(I , J)-R*CT(K ,J) 
CI(I,J)=CI(I,J)-R*CI(K,J):NEXT:NEXT:NEXT 
FOR K=1 TO N:D=CT(K,K) 

FOR 1=1 TO N:CT(K,I)=CT(K,I)/D:CI(K,I)=CI(K,I)/D:NEXT:NEXT 
PRINT: PRINT "THE RESIDUAL IN THE MATRIX THAT WAS INVERTED:" 

FOR 1=1 TO N:FOR J=1 TO N:PRINT USING "###. ######";CT(I , J) ; : NEXT: PRINT :NE 
PRINT: PRINT "INVERSE MATRIX:" 

FOR 1=1 TO N:FOR J=1 TO N:PRINT USING "##*.######”; Cl (I , J) ; : NEXT: PRINT :NE 
PRINT 

’7. COMPUTE THE TEMPERATURES FOR THE VARIOUS RIGHT-HAND SIDES 

'7.1. COMPUTE THE RHS 

'ON THE FIRST PASS, ALL HEATERS ARE OFF 

FOR 1=1 TO N:FOR K=1 TO M:NH(I ,K )=1 : NEXT: NEXT 

FOR K=1 TO M:F0R J=1 TO N:EP=0 :AL=0: IF A(J)=0 THEN 1230 

FOR 1=1 TO NC:EP=EP+AR(J,I)*EP(I):AL=AL+AR(J,I)*AL(I):NEXT 

R(J)=-QH(J ,K)-Q(J,K)-EP*A(J)*E(J,K)-AL*A(J)*S(J,K):NEXT 

'7.2. COMPUTE THE TEMPERATURES 

FOR 1=1 TO N:T(I ,K)=0:TU(I ,K)=0:FOR J=1 TO N:T(I ,K)=T(I,K)+CI(I, J)*R(J) 
NEXT: NEXT: NEXT: IF C$="2" THEN C$="3":GOTO 2170 'TO RECOMPUTE THE HTRS 
WO=W 

W=0:FOR K=1 TO M:R=0:FOR 1=1 TO N : R=R+QH ( I , K ) : NEXT : IF~ R W THEN W=R 
NEXT: FOR 1=1 TO N;W=W+.1*HC(I) :NEXT 
FOR K=1 TO M: PRINT "MISSION: ";K 

PRINT "NODE" ;TAB (6); "TEMPERATURE" ;TAB( 19); "HEATER PWR"; 

PRINT TAB (31); "HEATER CAP";TAB(43);"ALPH";TAB(49);"EP";TAB(54); 

PRINT "AREA RATIO BY COATING NO." 

PRINT TAB(6) ; "MIN (F) MAX";TAB(21);"BTU/HR ,, ;TAB(33);"BTU/HR";TAB(54); 
PRINT "1";TAB(60) ;"2";TAB(66) ;"3";TAB(72) ;"4" 

FOR 1=1 TO N: PRINT I;TAB(6); 

PRINT USING ”###.# ";T(I,K)-TU(I,K)-460;T(I,R)+TU(I,K)-460; 

PRINT TAB (19) ;QH(I,K) ;TAB(31 ) ;HC(I) ; 

IF A(I ) *0 THEN PRINT: GOTO 1440 

EP=0:AL=0:FOR J=1 TO NC:EP=EP+AR(I , J)*EP( J) :AL=AL+AR(I , J)*AL(J):NEXT 
PRINT TAB(43) ; :PRINT USING " .###" ;AL; : PRINT TAB(48); 

PRINT USING ".W;EP;:PRINT TAB(54); 

FOR J=1 TO NC: PRINT USING "#.*## ";AR(l, J) ; : NEXT: PRINT 
NEXT : PRI NT : NEXT K 

PRINT "PREVIOUS WEIGHT FUNCTION WAS" ;TAB(30) ;WO; 

PRINT TAB(40); "BTU/HR EQUIV AVERAGE HEATER POWER 
PRINT "THE NEW WEIGHT FUNCTION IS" ; TAB (30) ;W: PRINT 
'8. SELECT NEW SURFACE PROPERTIES OR NEW HEATER CAPACITIES 
IF L$="ON" THEN PRINT CHR$(15) 

PRINT "ENTER THE OPTION NUMBER TO BE EXERCISED:" 

PRINT " 1. ALTER THE HEATER CAPACITIES (COMPUTES REQ'D HEATER POWERS) 
PRINT " 2. ALTER THE SURFACE PROPERTIES (RESETS HEATERS TO ZERO) 

PRINT " 3. ALTER THERMOSTAT SET POINTS (COMPUTES REQ’D HEATER POWERS) 

PRINT " 4. COMPUTE THE DERIVATIVES 

PRINT " 5. TOGGLE LINE PRINTER 


A-9 



TABLE A-3 (Continued) 

THERMAL OPTIMIZATION PILOT PROGRAM 


01560 PRINT " 6. STORE THE INPUT DATA ON DISK 

01570 PRINT " 7. LIST THE TEMPERATURE DATA AGAIN 

01580 PRINT " 8. INCLUDE THE TOLERANCES 

01590 PRINT " 9. END 

01600 INPUT "WHICH";C$:IF C$="2" THEN 1970 

01610 IF C$="7" THEN 1300 

01620 IF C$="9" THEN END 

01630 IF L$="0N" THEN PRINT CHR$(14) 

01640 IF C$ "6" THEN 1770 

01650 PRINT "THE DATA WILL BE STORED IN FILE NAME ";S$;:INPUT " OK (Y/N)";C$ 
01660 IF C$="Y"THEN 1680ELSE INPUT "DO YOU WANT TO RETURN TO THE MENU (Y/N)";C$ 
01670 IF C$="N" THEN INPUT "WHAT SHOULD THE FILE NAME BE";S$ ELSE 1490 

01680 IF S$="N" THEN PRINT "CAN'T USE N FOR A FILE NAME": GOTO 1490 

01690 OPEN "0",1,S$ 

01700 PRINT#1,N,M,NC 

01710 FOR 1=1 TO NC:PRINT#1 ,NC$(I) ,AL(I) ,EP(I) :NEXT 
01720 FOR 1=1 TO N-l:FOR J=I+1 TO N:PRINT#1 ,C(I , J) :NEXT:NEXT 
01730 FOR 1=1 TO N:PRINT#1 ,A(I) ,T(I , 1 ) :IF A(I)=0 THEN 1750 
01740 FOR J=1 TO NC: PRINT# 1 ,AR(I , J) :NEXT 

01750 NEXT : FOR 1=1 TO M:FOR J=1 TO N: PRINT# 1 ,Q(J,I),S(J,I),E(J,I) : NEXT: NEXT 

01760 FOR 1=1 TO N:PRINT#1 ,SP(I) : NEXT: CLOSE: GOTO 1490 

01770 IF L$="OFF" AND C$="5" THEN L$="ON":PRINT CHR$(14):GOTO 1490 

01780 IF L$="ON" AND C$="5" THEN L$="OFF": PRINT CHR$ ( 1 5 ) : GOTO 1490 

01790 IF C$="3" THEN 2090 

01800 IF C$="8" THEN 2575 

01810 IF C$="4" THEN 2760 

01820 IF C$ "1" THEN PRINT "OPTION NOT YET AVAILABLE": GOTO 1490 

01830 '8.1. SELECT THE HEATER CAPACITIES 

01840 PRINT: PRINT "ENTER THE HEATER CAPACITIES:" 

01850 IH=0 

01860 PRINT "NODE"; TAB (10); "SET POINT (F)";TAB(25) ;"BTU/HR CAPACITY" 

01870 C$="N":FOR 1=1 TO N: PRINT I ; TAB ( 10) ;SP(I) ;TAB(25) ;HC(I) ;TAB(35) ;CH$; 

01880 R=HC(I ) : INPUT HC(I):IF HC(I)= R THEN HC( I )=0: IH=IH+1 : IH(IH)=I:GOTO 1910 

01890 FOR K=1 TO M:NH(I,K)=2:IF QH(I,K) HC(l) THEN QH(I ,K)=HC(l) 

01900 NEXT 
01910 NEXT 

01920 FOR K=1 TO M:FOR J=1 TO N:IF A(J)=0 THEN EP=0:AL=0:GOTO 1940 
01930 EP=0:AL=0:FOR 1=1 TO NC:EP=EP+AR(J,I)*EP(I):AL=AL+AR(J,I)*AL(I):NEXT 

01940 FOR J=1 TO N:R(J)=-QH(J,K)-Q(J,K)-EP*A(J)*E(J,K)-AL*A(J)*S(J,K):NEXT 

01950 FOR 1=1 TO N:T(I,K)=0:FOR J=1 TO N:T(I,K)=T(I,K)+CI(I,J)*R(J):NEXT:NEXT 
01960 NEXT: GOTO 2220 

01970 '8.2. SELECT THE NEW SURFACE PROPERTIES 

01980 PRINT "NODE" ;TAB (6) ; "ALPHA"; TAB (13); "EPSILON" ;TAB (25); "AREA FRACTION FOR" 
01990 PRINT " COATINGS": PRINT TAB(25); 

02000 FOR J=1 TO NC:PRINT USING ".### ";AL(J); :NEXT:PRINT "(ALPHA)": PRINT TAB 

02010 FOR J=1 TO NC: PRINT USING ".### ";EP(J); :NEXT:PRINT "(EPSILON)" 

02020 FOR 1=1 TO N:IF A(l)=0 THEN 2070 

02030 PRINT I;TAB(6) ;AL;TAB(13) ;EP;TAB(25) ; :FOR J=1 TO NC 
02040 PRINT USING "#.### ";AR(I,J); :NEXT:PRINT " CHANGE (Y/N)"; 

02050 X$="N": INPUT X$:IF X$ "Y" THEN 2080 

02060 PRINT TAB (25); : INPUT AR(I , 1 ) ,AR(I ,2) ,AR(I,3) ,AR(I ,4) 

02070 HC(I)=0:F0R K=1 TO M:QH(I ,K)=0:NEXT 

02080 NEXT : GOTO 840'BECAUSE THE DIAGONAL (VIA EP) IS CHANGED 
02090 '8.3. COMPUTE THE REQUIRED HEATER CAPACITIES 
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02100 

02110 

02120 

02130 

02140 

02150 

02155 

02160 

02170 

02180 

02190 

02200 

02210 

02220 

02230 

02240 

02250 

02260 

02270 

02280 

02290 

02300 

02310 

02320 

02330 

02340 

02350 

02360 

02370 

02380 

02390 

02400 

02410 

02420 

02430 

02440 

02450 

02460 

02470 

02480 

02490 

02500 

02510 

02520 

02530 

02540 

02550 

02560 

02570 

02575 

02580 

02590 

02600 

02610 


'BECAUSE THE TEMP UNCERTAINTIES ARE NOT LINEAR BUT RSS, THE HEATER WILL 
'NOT MAKE THE TEMPS EXACTLY MEET THE THERMOSTAT SET POINT, AT LEAST AS 
’WE NOW SOLVE THE PROBLEM. WE TREAT THE UNCERTAINTY AS A LINEAR EFFECT. 
PRINT "ENTER THE THERMOSTAT SET POINTS : M sIH=0: FOR 1=1 TO N:HC(l)=0:NEXT 
PRINT "NODE" ; TAB ( 10) ; "SET POINT (F) (SP=-999 = NO THERMOSTAT OR HEATER) 

FOR 1=1 TO N:SP(l)=SP(l )-460:PRINT I;TAB(10) ;SP(l) ;TAB(25) ;CH$; :INPUTSP(l 
SP(I )=SP(I)+460:NEXT 
PRINT 

FOR 1=1 TO N:IF SP(I) -999 THEN FOR K=1 TO M:NH(I,K)=0:NEXT 
NEXT 

FOR 1=1 TO N:HC(I)=0:NEXT 

'TO INCORPORATE THE EFFECTS OF EACH HEATER ON ALL OF THE TEMPERATURES, 

'WE MUST SOLVE A MATRIX EQUATION FOR EACH NEW HEATER POWER (BY MISSION) 
FOR K=1 TO M 

IH=0:F0R 1=1 TO N:IF NH(l,K) 0 THEN 2250 

IH=IH+ 1 : IH ( IH ) =1 

NEXT 

IF IH=0 THEN 2440 

IF IH=1 THEN I=IH(1) :CN(1 , 1 )=1/CI(I ,1) :GOTO 2370 

FOR J=1 TO IH : JH=IH ( J ) : FOR L=1 TO IH:LH=IH(L) :CT( J,L)=CI (JH ,LH) 

CN( J,L)=0: NEXT: CN( J ,J)=1 :NEXT 

FOR KK=1 TO IH-1 :FOR II=KK+1 TO IH:R=CT(II ,KK)/CT(KK ,KK) :FOR JJ=1 TO IH 
CT(II , JJ)=CT(II , JJ )-R*CT(KK, JJ) :CN(II,JJ)=CN(II , JJ)-R*CN(KK , JJ ) : NEXT: NEXT 
NEXTiFOR KK=IH TO 2 STEP -l:FOR II=KK-1 TO 1 STEP -1:R=CT(II,KK)/CT(KK,KK 
FOR JJ=IH TO 1 STEP -1 : CT(II , JJ)=CT(II , JJ)-R*CT(KK , JJ) 

CN(II , JJ )=CN(II , JJ )-R*CN(KK , JJ ): NEXT : NEXT : NEXT 

FOR KK=1 TO IH:D=CT(KK,KK):FOR 11=1 TO IH 

CT(KK ,11 )=CT(KK , II )/D: CN(KK , II )=CN(KK , II )/D: NEXT: NEXT 

FOR 11=1 TO IH: I=IH(II) :R(II)=0:FOR J=1 TO IH:JH=IH(J) 

R(II)=R(II)-CN(II,J)*(SP(JH)-T(JH,K)+TU(JH,K)):NEXT 

IF QH(I,K)+R(II)= 0 THEN 2410 

QH(I,K)=0:R(II)=0:NH(I,K)=1:GOTO 2230 

NEXTiFOR 11=1 TO IH: I=IH(lI);QH(l,K)=QH(I,K)+R(II) 

IF QH ( I , K ) HC(l) THEN HC(I)=QH(I ,K) 

NEXT 
NEXT K 

GOTO 1210 'BECAUSE T'S NEED TO BE RECALCULATED BEFORE ANY OTHER CHANGES 
'WE NEED TO EXAMINE HOW TO INCORPORATE TOLERANCES 

'PROBABLY WE SHOULD USE THE RMS ERROR AND GIVE THE TEMP OR HTR RANGE 
'THE USER COULD SPECIFY THE NUMBER OF SIGMAS TO CONSIDER, WITH THE 
'INPUT BEING ENTERED AS 3 SIGMA. THE SET POINTS ENTERED BY THE USER 
'SHOULD BE TREATED AS THE AVG - N*SIGMA WHEN COMPUTING THE HEATER 
'POWER. IN PRACTICE, WE MAY NOT ALWAYS BE ABLE TO PUT THE HEATER ON 
'THE CRITICAL NODE. WE HAVE NOT ALLOWED FOR THAT IN THE OPTIMIZATION. 

'WE SHOULD PROBABLY ADD FIXED-TEMP NODES TO THE PROBLEM STATEMENT, 
'ALTHOUGH IT IS DIFFICULT TO CONCEIVE OF THE PHYSICAL SITUATION. 

'8. INCLUDE THE TOLERANCES 

'FOR NOW, WE ASSUME THAT THE COATING PROPERTIES ARE KNOWN TO WITHIN 
*0.03, REGARDLESS OF THE COATING. THE TOLERANCES ARE DESIGNATED AU,EU 
PRINT:C$="IF COATING TOLERANCES ARE INCLUDED": PRINT TAB(40-LEN(C$ )/2) :C$ 
FOR J-l TO NC: AU( J)=.03:EU( J )=.03:NEXT 

FOR K=1 TO M:FOR J=1 TO N:TU(J,K)=0:NEXT:FOR JC=1 TO NC:L=0:FOR J=1 TO N 
IF A(J)=0 OR AR( J , JC)=0 THEN R(J)=0:GOTO 2630 
R(J)=AR(J,JC)*A(J)*(E(J,K)-H(J)*T(J,K))*EU(JC) 
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02620 L-l 

02030 NEXT J 

02640 IF L-0 THEN 2740 

02650 'SOLVE FOR ALL OF. THE TEMPERATURE UNCERTAINTIES FIRST, AS IF NONLINEAR 

02660 FOR J=1 TO N:CT(J,1)-0:F0R 1-1 TO N:CT(J,1)-CT(J,1)+CI(J,I)*R(I):NEXT:NEX 

02670 FOR J-l TO N:TU(J,K)«TU(J,K)+CT(J,1)*CT(J,1):NEXT 

02680 'NOW INCLUDE THE EFFECTS OF ALPHA 

02690 FOR J=1 TO N:IF A(J)-0 OR AR(J,JC)-0 THEN R(J)»0:GOTO 2710 
02700 R(J)-AR(J,JC)*A(J)*S(J,K)*AU(JC) 

02710 NEXT 

02720 FOR J-l TO N:CT(J,1)«0:FOR 1-1 TO N:CT(J,1)-CT(J,1)+CI(J,I)*R(I):NEXT:NEX 

02730 FOR J-l TO N:TU(J,R)-TU(J,K)+CT(J,1)*CT(J,1):NEXT 

02740 NEXT JCsFOR J-l TO N:TU(J,R)-SQR(TU(J,R)):NEXT:NEXT RsGOTO 1300 

02750 'COMPUTE THE DERVIATIVES WITH RESPECT TO THE AREA RATIOS 

02760 FOR K-l TO M ’DE(J,K,I,L) - D(T(J,K))/D(AR(I,L)) 

02770 FOR 1-1 TO N:IF A(I)»0 THEN 2900 

02780 FOR L-l TO NC 

02790 R(I)=A(I)*(AL(L)*S(I,R)+EP(L)*(E(I,R)-H(I)*T(I,R))) 

02800 FOR J-l TO N:DE(J,R,I,L)»CI(J,I)*R(I):NEXT 

02810 NEXT L 

02820 C$«"MISSION"+STR$(R) 

02830 PRINT: PRINT TAB(40-LEN(C$)/2);C$ 

02840 PRINT "DERVIATIVES OF ALL NODE TEMPERATURES WITH RESPECT TO THE AREA RATI 
02850 PRINT "ON N0DE"jI;"F0R COATING:";TAB(27);"l";TAB(37);"2";TAB(47);"3"; 
02860 PRINT TAB (57); "4" 

02870 PRINT TAB (10); "AFFECTED NODE";TAB(32) ;"D ERIVATIVES" 

02880 FOR J-l TO N:PRINTTAB( 15) ; Jj :F0R L-l TO NCsPRINT TAB(i7+L*10);DE(J,R,I,L) 
02890 NEXT SPRINT: NEXT 

02900 NEXT I SPRINT: NEXT K 

02910 GOTO 1490 
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